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We study a symmetry property of momentum distribution functions in the steady state 
of heat conduction. When the equation of motion is symmetric under change of signs for 
all dynamical variables, the distribution function is also symmetric. This symmetry can be 
broken by introduction of an asymmetric term in the interaction potential or the on-site 
potential, or employing the thermal walls as heat reservoirs. We numerically find differences 
of behavior of the models with and without the on-site potential. 
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1. Introduction 

Recently, considerable progress has been made in statistical mechanics of nonequilibrium 
steady states beyond the linear-response regime. First, new formulas such as the fluctuation 
theorem 1 and the Jarzynski equality 2 have been found. These formulas can be used as a 
basis for arguments on transport phenomena in small systems, such as energy transfer in mo- 
tor proteins or heat and electric transport in nano-scale devices. Second, Sekimoto proposed 
energetics on Langevin systems, namely systems in the scale where fluctuations become im- 
portant. 3 It was further developed and some interesting results were obtained 4 Lastly, some 
exactly solvable models were found, 5 ' 6 which give us important insights to the problem. 

Lattice thermal conduction is a typical example of the nonequilibrium steady states. It 
may be more fundamental than stochastic systems because it is described by a Hamiltonian. 
Study of thermal conductivity has a long history, but many important results were obtained 
only recently. It was established already in 1960s that integrable systems do not have tem- 
perature gradient and heat flux is proportional to not temperature gradient but temperature 
difference between the ends. For example, heat conduction in the harmonic chain was solved 
in 1958. 7 However, it was found in 1997 that the thermal conductivity in the Fermi-Pasta- 
Ulam (FPU) chain diverges like k ~ N a (a ~ 0.37) as the system size N is increased. 8 Since 
then, a lot of researches were carried out and it was clarified that the divergence is caused 
by a long-time tail in current autocorrelation function, which is originated from momentum 
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conservation. 9 ' 11 ' 24 On the contrary, if the system includes a potential force from substrate 
the thermal conductivity is convergent. The (j) 4 system, 10 the ding-a-ling model 12 and the 
ding-dong model 13 belong to the class of systems. However, we must mention that exceptions 
exist for this rule and that proposed theories have some discrepancy concerning the value of 
exponent a. 11,14 Thus, we are still far from complete understanding. 

As we see in the above, the study of lattice thermal conduction has been restricted to 
the Fourier heat law and existence or nonexistence of thermal conductivity. In fact, there 
are few studies on distribution functions that describe the nonequilibrium steady state of 
heat conduction. To construct nonequilibrium statistical mechanics, however, properties of 
the distribution function should be studied. Thus, in this paper we investigate characteristics 
of the nonequilibrium steady states. 

Analogously with the Gibbs ensemble in equilibrium statistical mechanics, there should 
be a distribution function on the phase space that describes a nonequilibrium steady state. 
However, the phase space is high-dimensional and numerically intractable. So, we focus on 
momentum distribution of a single particle. In equilibrium, it is a Maxwellian distribution at 
some temperature. We study how it deviates from the Maxwellian distribution in nonequlib- 
rium steady states. 

In particular, we focus on a symmetry property. In equilibrium systems, the momentum 
distribution does not depend on the potentials. On the other hand, symmetry of the potentials 
can affect the momentum distribution in nonequilibrium states. We numerically investigate 
what happens if the symmetry is broken by introducing a small asymmetric term into the on- 
site potential or the interaction potential or employing the thermal wall as the heat reservoir. 
As a result, differences are found in the behavior of the models with and without on-site 
potentials. This may be relevant with the reported behavior of heat conduction and might be 
useful for full understanding of the problem. 

In fact, Aoki and Kusnezov 16 discussed similar deviation from the Maxwellian distribution 
and derived a scaling form for the fourth-order cumulants. However, they were limited to 
symmetric systems and symmetric deviations. In this paer, we extend to asymmetric models 
and asymmetric deviations. 

In Section 2, we describe the models and heat reservoirs employed in our simulations. In 
Section 3, we demonstrate a relation between the symmetry in the equation of motion and that 
in the momentum distribution functions. In Section 4, we show our numerical results when a 
weak asymmetry is introduced to a symmetric model. Section 5 is devoted to summary and 
discussion. 
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We consider one-dimensional systems composed of N particles of unit mass with a Hamil- 
tonian of the form 



N 



n=l 



H({q n },{ Pn }) = J2 ^ + U(q n ) + 5>(< ?n+1 -<?„), (1) 



2 



N 



n=0 



where q n and p n are the displacement and the momentum of particle n £ {1,2, ... , N}, U(q) is 
an on-site potential, and V(q) is an interaction potential between nearest-neighbor particles. 
We impose the fixed boundary condition, which is represented by setting go = qN+i = in Eq. 
(1). Various models are generated by varying U and V. Typical examples are (a) the Fermi- 
Pasta-Ulam (FPU) model (U{x) = 0, and V(x) = ^ + ^), (b) the (j) A model (U(x) = 4 + 

2 

and V(x) = ^-), and (c) the Toda model (U(x) = and V(x) = exp(x) — x). 

The particles at the ends of the chain are in contact with two heat reservoirs: one at 
temperature Tl on particle 1 and the other at temperature Tr on particle N. If Tl = Tr, the 
system goes to equilibrium, and if Tl / Tr, heat conduction occurs. Thus, for 2 < n < N — 1, 
the equations of motions are given as 

q n =Pn, Pn = -U'(q n )+V'(q n+1 -q n )-V'(q n -q n - 1 ), (2) 

while those for particles 1 and N are modified from the Hamiltonian form. In this paper, 
we consider the following three kinds of reservoirs to study what differences are generated 
by the types of heat reservoirs. The first is the Langevin reservoir that is given by adding a 
dissipation term and fluctuating force to the equation of the motion as 

Pi = ~U'( qi ) + V'(- gi ) - V\ gi - q 2 ) ~ IPX + (3) 

and 

PN = ~U'(qN) + V'(qN-i ~ qN) ~ V(qN) - JPN + (4) 

In the above equations, and £r(*) are Gaussian white noise with correlation 

(UtW)) = 2ik B T a 8 aP 8{t - t') 

where a, (3 = L or R. The second is the Nose-Hoover reservoir, 26 which is a kind of determin- 
istic reservoir and widely used in the literature. In this reservoir, the equations of motion for 
particles 1 and iV are modified as 

Pi = -U'(q 1 ) + V'(q 2 -q 1 )-V'(q 1 )-t; L p 1 , & = 1 ^ - l) , (5) 



and 



PN = -U'{q N )+V'{-q N )-V'{q N -q N _ 1 )-ti R p N , ( R = ^ (j| - , (6) 
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and time evolution of £l and £r are given by 
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and 



l 
e 
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(8) 



where, 9 is a constant which means response time of the thermostat. The last is the thermal 
wall. 15 When a particle hit the thermal wall, it is reflected with a momentum chosen according 
to the probability distribution 



where T denotes the temperature of the wall. For our systems, we place a thermal wall at 
temperature Tl on the left hand of particle 1 and the other wall at temperature Tr on the 
right hand of particle N . 

The three models mentioned in Section 1 are representatives of the types of heat conduc- 
tion behavior. Namely, the FPU model represents the models where the thermal conductivity 
diverges in the thermodynamic limit due to a long-time tail in the autocorrelation function 
of heat flux. The (j) 4 model represents the class of models where the thermal conductivity 
converges as the system size is increased. The Toda model represents the integrable models 
where temperature gradients are not formed and heat flux is proportional to temperature 
differences between the reservoirs. 

3. Deviation from the Maxwellian distribution 

We carriy out numerical simulations on the three models with Langevin heat reservoirs 
at Tl = 2.0 and Tr = 1.0 and calculate the single-particle momentum distribution functions 
Pn(p) {n = 1, . . . , N) as follows. The deterministic part of the equation of motion is solved 
with a mixed use of the sixth-order symplectic method and the fourth-order Runge-Kutta 
method with a single time step At = 0.01. After computing the deterministic part, we add a 
Gaussian random noise with variance 2^Tl At to the momentum of particle 1, and one with 
variance 2jTjiAt to the momentum of particle TV. To compute the distribution functions, 
simulation was done during 10 10 -steps after the system reached a steady state. Figure 1 shows 
the difference between the computed momentum distribution function and the Maxwellian 
distribution function 



where temperature T n is determined via average of kinetic energy, namely T n = 

We find that the deviation is symmetric in the FPU model and the model, while it is 
asymmetric in the Toda chain. This is because the equations of motion are symmetric under 
the transformation of changing signs of the variables for the former two models, but not for 




(9) 




(10) 
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Fig. 1. Deviation from the Maxwellian distribution for (a) the FPU model, (b) the (f) 4 model, and (c) 
the Toda model. Here the system size is N = 32 and the temperature values of the heat reservoirs 
are given as Tl = 2.0 and Tr — 1.0. 10 10 iterations are carried out for (a), (b) and (c). 



the Toda chain. In general, the equation of motion (2) is invariant under the transformation of 
changing the signs of the variables . . . ; C[n->Pn) — ► {—Qi-, —Pi] —Qn, ~Pn) if both the 

on-site potential U(q) and the interaction potential V(q) are even functions, i.e., U(—q) = U(q) 
and V(— q) = V(q). We here note that there is no physical reason that the interaction potential 
must be even. Denoting lattice constant by c, equality V(—q) = V(q) means that two particles 
with distances c + q and those with distance c — q have the same interaction energy, which is 
not expected in general. 

Moreover, the Langevin equations for particles 1 and N, (3) and (4), are also in- 
variant under the transformation if the signs of the noise terms and £r(<) are 
also changed. The last operation does not change the statistics of the noise terms. Simi- 
larly, Nose-Hoover reservoir is also symmetric under operation (qi,p±; . . . ; qN,PN] £l, £r) — ► 
(— <Zi, — Pi; ■ ■ ■ ; — 9jv, — Pjv;£i,)£r)- On the other hand, since the thermal wall introduces an 
asymmetric potential, it breaks the symmetry. 

The symmetric deviation seen in the FPU model and the 4> 4 model is in contrast with the 
antisymmetric deviation seen in particle systems described by the Boltzmann equation. 27 In 
the latter case, the expectation value of heat current is represented as \ J pp 2 f(p)dp with use 
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of the momentum distribution function f(p). Clearly, it vanishes if f(p) is symmetric. Thus, 
antisymmetric deviation must exist. On the other hand, in our lattice models the expectation 
value of heat current is represented as 

- J (Pn+Pn+i) V'(q n+ i - q n )Pn,n+i{q n ,Pn, q n +i,Pn+i)dq n dp n dq n+1 dp n+1 , 

where P n ,n+i{Qn,Pm Qn+i,Pn+i) denotes a two-body distribution function in the steady 
state. Because V'(q) is antisymmetric when V(q) is symmetric, the integral does not 
vanish if P n ,n+i(qn,Pn, q n+ i,p n+ ±) is symmetric. This also indicates the existence of 
correlation. Namely, if P n ,n+i(<ln,Pn,qn+i,Pn+i) is decomposed into a product form 
Pn,n+i{q n ,Pn,qn+i,Pn+i) = P^l+i (<?n, Qn+i)P^l +1 {Pn, Pn+i) as in the equilibrium case and 
Pnn+i and/or P^\ + i is symmetric, the heat current vanishes. Thus, to produce nonzero heat 
current some correlation must exist between momentum and position. The symmetric devia- 
tion in the momentum distribution function is a reflection of such correlation. 

4. Properties of asymmetric deviation 

In this section, we examine how the momentum distribution functions are affected if 
the symmetry is broken by small modifications to a symmetric system. As we noted in the 
Introduction, there are three kinds of such modifications, which we describe in the following. 
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Fig. 2. The third-order moments (p^) in the system with asymmetric interaction potential (11) and 
on-site potential U(q) = 0. The system size is N — 16, parameter values are a = —0.3, 0, 0.3, 0.5, 
0.7, 1.0. Langevin reservoirs at Tl = 2.0 and Tr = 1.0 are used. 

The first is the case of modifying the interaction potential. Let us consider a system with 
the following interaction potential, 

V(q) = \q 2 + ^ 3 + \q\ (11) 

with a constant a representing the magnitude of asymmetry. The on-site potential U(q) is set 
to be zero. We carry out numerical simulations of the system of size N = 16 and calculate 
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the third-order moments 

/+oo 
p 3 P n (p)dp (n = l,...,N). (12) 
-oo 

The results are shown in Fig. 2, which indicates that asymmetric deviations from the 
Maxwellian distribution emerge in the whole system. Namely, the moments are in the same 
order of magnitude except a few particles near the ends. This behavior is not changed if the 
system size is varied. See Fig. 3. Although there appear relatively large variations near both 
the ends, changes are smooth in the middle of the system. Figure 4 shows that the third- 
order moment changes in proportion to a when a is small. But the tendency changes around 
a ~ 2, where the moment have an extremum. This is because the stability of q = is lost at 
a = 2 and two stable points with q ^ appear beyond a = 2. Then, q loses the meaning of 
displacement from the stable point in the latter case. 
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Fig. 3. The third-order moments | (p^) | in various sizes of the system having the asymmetric potential 
(11) with a = 1. The reservoir temperatures are given as Tl = 2.0 and Tr = 1.0. 
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Fig. 4. Variation of {p\) as the parameter a is varied. The system is the same as in Fig. 2 



We check this behavior with some conditions modified. One is the introduction of on-site 
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potential. When we use the nonlinear on-site potential 

u( q ) = \e + \q\ 

the third-order moments become almost constant except for the end particles and this behavior 
does not change with the system size. The next is the use of Nose-Hoover reservoirs instead of 
Langevin reservoirs. In this case, the magnitude of the third-order moments is larger than the 
Langcvin case and their spatial variation is rather smooth. No noticeable peaks are formed 
if the system size is changed. The next is the variation of the boundary temperatures. When 
the temperatures of the reservoirs are given as Tl = 60 and Tr = 30, the overall behavior is 
the same as the low temperature case, only magnifying the absolute values of the moments. 
The last is insert of the fifth-order term instead of the third-order term into the interaction 
potential. The use of the fifth-order term leads to little difference from the case of third-order 
term. In all cases, the overall appearance of the third-order moments are observed. 
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Fig. 5. Deviations from the Maxwellian distribution for (a) particle 8 and (b) particle 16 in the 
model with the asymmetric on-site potential U(q) — \q 2 + ^q 3 + jq 4 and the harmonic interaction 
potential V(q) = \q 2 . The system size is N = 16 and Langevin reservoirs at Tl = 2.0 and Tr = 1.0 
are used. 10 10 iterations are carried out for each of them. 



Next, we consider the case of modifying the on-site potential. Then, the third-order term 
is inserted into the on-site potential of the <p A model as 

U( q ) = \f + If + \d\ (13) 

and the harmonic interaction potential V(q) = q 2 /2 is used. Figures 5 (a) and (b) show the 
deviation from the Maxwellian distribution for particles 8 and 16 in the system with a = 1.0 
of system size N = 16. Asymmetric deviation is observed for the end particle but invisible in 
the middle of the chain. Figure 6 shows the third-order moments (p^) and its variation with 
parameter a, which clearly indicates that only the particles at the ends of the chain are affected 
by the asymmetry. Notice that even for the end particles the magnitude of moments is much 
smaller than that in the case of modifying interaction potential. Namely, asymmetry in the 
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Fig. 6. The third-order moments (p^) in the system with the asymmetric on-site potential U(q) = 
^- + a^- + 2j- and the harmonic interaction potential V(q) = The system size is N = 16, 
Langevin reservoirs at Tl = 2.0 and Tr = 1.0 are used, and parameter values are a = —0.3, 0, 
0.3, 0.5, 0.7, and 1.0. 
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Fig. 7. Variation of (pf), (p%/ 2 ) and (p^) with parameter a. The system and reservoir are the same 
as in Fig. 6 

on-site potential produces much smaller effects than asymmetry in the interaction potential. 
Figure 7 shows (p^ t ) vs a for particles at the ends and one in the bulk. The moments change 
linearly with a for particles at the ends but shows little changes for the particle in the bulk. 
The linearity implies that this change is certainly caused by the introduced asymmetry. In 
the middle of the system, the effect of asymmetry is too small to be observed. 

When the interaction potential is not a harmonic potential but a nonlinear one as 

V(q) = \q 2 + \q\ (14) 

the largest deviation of the momentum distribution is observed not for particles 1 and N but 
particles 2 and N—l. However, the effect of asymmetry is still confined to the vicinity of both 
the ends. Figure 8 (a) shows the third-order moments near the right end and their system-size 
dependence. As seen from this figure, the number of particles affected by the asymmetry is 
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not changed with the system size N. If we further replace the Langevin reservoirs with the 
Nose-Hoover ones, the confinement of the asymmetric effects is not changed, as is shown in 
Fig. 8 (b). A notable change is that the magnitude of the moments becomes fairly large. The 
shape of the deviation depends on the interaction potential and the reservoirs. Figure 9 (a) 
shows the deviation of the momentum distribution for particle N — 1, which has the largest 
moment in the system with the nonlinear interaction potential and Langevin reservoirs. The 
deviation has similar shape with Fig. 5 (b). However, the deviation changes its sign if the 
reservoirs are replaced with the Nose-Hoover ones as seen in Figure 9 (b), which represents 
the deviation for the rightmost particle. In this case, the deviation for particle N — 1 and that 
for particle N have different signs. From these detailed observations, we conclude that the 
asymmetry effect is confined to the vicinity of the ends in the systems with an asymmetric 
on-site potential and a symmetric interaction potential. On the other hand, the extent of the 
confinement and the shape of the deviation depend on the details of the systems. 




Fig. 8. The third-order moments (p^) near the right end in the system with the asymmetric on-site 

2 3 4 2 4 

potential U(q) = \ + ^- + \ and the interaction potential V(q) = \ + \- The system size 
N = 16, 32, 64, 128, 256 and 512. Heat reservoirs are (a) Langevin reservoirs and (b) Nose-Hoover 
reservoirs at Tl = 2.0 and Tr = 1.0. 



Now we turn to the case where the heat reservoirs break the symmetry. As mentioned in 
Sec. 3, the Langevin reservoirs and the Nose-Hoover reservoirs keep the symmetry but the 
thermal wall does not. Figure 10 compares the profiles of the third-order moments in (a) the 
FPU model and (b) the <j) A model when the thermal walls are employed as heat reservoirs. 
In the (p A model, the influence of the symmetry-breaking is virtually limited to only the four 
particles, 1, 2, N — 1 and N. The moments almost vanish for the other particles. On the 
contrary, in the FPU model, finite moments appear even in the central part of the system, 
though the magnitude is much smaller than the case of particles near the boundaries. Figures 
11 (a) and (b) enlarge the profile of moments near the ends, and show system-size dependence 
for the FPU model and the 4> 4 model, respectively. In the <p A model, the number of sites 
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Fig. 9. Deviations from the Maxwellian distribution in the system of size 16 with the asymmetric 
on-site potential U(q) = \q 2 + 1<? 3 + jq 4 and the nonlinear interaction potential V(q) = \q 2 + \q A - 
(a) Particle 15 in the system with Langevin reservoirs, (b) Particle 16 in the system with Nose- 
Hoover reservoirs. The reservoir temperatures are given as Tl = 2.0 and Tr = 1.0. 10 10 iterations 
are carried out for each of them. 
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Fig. 10. The third-order moments (p 3 ) for (a) the FPU model and (b) the <f> model when the thermal 
walls at Tl = 2.0 and Tr = 1.0 are used. The system size is iV = 32. 



affected by the asymmetry does not change. In the FPU model, the moments behave as if 
to decrease exponentially from the boundary into the bulk. Interestingly, however, it actually 
goes beyond zero and finite moments are obtained in the central region. Recall that the (j) 4 
system has a bulk limit of the thermal conductivity, while the FPU model does not. Our result 
also implies that the FPU system has no bulk. 

In the dilute gas system between two thermal walls at different temperatures, it is known 
that asymmetric deviations appear in the whole system. 2 Thus the lattice system is totally 
different from the gas system in this point. 

5. Summary and Discussion 

We have numerically studied the single-particle momentum distribution functions in one- 
dimensional lattice dynamical system in nonequilibrium steady states of heat conduction. We 
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Fig. 11. The third-order moments (p„) for (a) the FPU model and (b) the (j) 4 model when the thermal 
walls at T L = 2.0 and T R = 1.0 are used. The system size is N = 16, 32, 64, 128, 256, 512. 



especially focus on the deviations from the Maxwellian distribution. This deviation reflects a 
symmetry of the system. Namely, if the system is invariant under the change of signs of the 
dynamical variables, the momentum distribution function must be even. This symmetry can 
be broken by the interaction potential or the on-site potential or the heat reservoir. In the first 
case, the effect of asymmetry extends to the whole system. In the second case, the effect of 
asymmetry is concentrated near the ends. In the last case, we have found notable differences in 
systems with and without on-site potentials. Namely, the effect of asymmetry is localized near 
the ends in the systems with on-site potentials, whereas the effect is extended to the center of 
the system without on-site potentials. This means there is no bulk limit in the system without 
on-site potentials, which is consistent with the known results on the convergence or divergence 
of thermal conductivity. 

We need to develop a theoretical explanation for our present results and clarify a relation 
between the deviation and thermal conductivity. Before studying that, it is necessary to specify 
what parameters are relevant to the deviations. As we mentioned in Sec. 3, deviation of the 
momentum distribution is related to two-site correlations including heat flux. Thus, heat flux 
and local temperature, which is the second-order moment of the distribution, are considered 
to be important parameters. It is a future problem to clarify if they are sufficient to determine 
the deviations. 
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